by Paulo Marques, 2014/03/09 (Adapted in 2018/10/15 to Python from Julia)
Since high school I've been fascinated with $\pi$ -- this infinite non-repeating irrational transcendent number. In fact, not only was I fascinated with $\pi$ but I was fascinated about everything related to it. In 11th grade I asked my math teacher about how to deduce the area of a circle. Her answer: for that you need to learn how to integrate – wait until university. But I couldn't wait and head to the library where I found a single book that talked about the subject – an obscure calculus book by an author named Piskunov. And, to integrate I've learned – just because of $\pi$. But I digress. ..
This story is not about calculus or "symbolic" integration. It's about how Archimedes calculated $\pi$ circa 200 B.C. In the "Measurement of a Circle" Archimedes states that:
"The ratio of the circumference of any circle to its diameter is greater than $3\tfrac{10}{71}$ but less than $3\tfrac{1}{7}$"
This is the first really accurate estimation of $\pi$. I.e., he calcuated $3.140845070422535 < \pi < 3.142857142857143$. A good approximation of $\pi$ is 3.141592653589793. So, this is two decimal places correct. That's pretty impressive.
According to the story, Archimedes did this by inscribing and circumscribing a circle with regular polygons and measuring their perimeter. As the number of sides increases, the better these polygons approximate a circle. In the end Archimedes was using a 96-sided polygon. The next image illustrates the idea.
One of the annoying things when books talk about this is that they always show this nice picture but never ever do the actual calculation. So, using Python how can we do this?
Let's start by assuming that we are going to use a circle with a radius of 1 and we inscribe a square in it. (The square's side is going to be $\sqrt{2}$.)
Now, assume that the side of this polygon is $s_n$ and you are going to double the number of sides where the length of each new side is $s_{n+1}$. We can draw several triangles in the figure that will help us out:
If we take the side $\overline{AB}$, which measures $s_n$, and break it in two, we get the triangle $\overline{ACD}$. This triangle has a hypotenuse of $s_{n+1}$, an adjacent side of $s_n/2$ and a height of $h$. Note that the new polygon that we are forming is going to have eight sides (i.e., double the number of sides we had), each one measuring $s_{n+1}$. From this we can write:
$$ h^2 + (\frac{s_n}{2})^2 = s_{n+1}^2 $$Looking at the triangle $\overline{BCO}$, which is rectangle, we note that: its hypotenuse is 1, one side measures $1-h$ and the other measures $s_n/2$. Thus, we can write:
$$ (1-h)^2 + (\frac{s_n}{2})^2 = 1^2 $$These two relations will always apply as we contantly break the polygons into smaller polygons. As we progress, the perimeter of the polygon $P_n$, obtained after $n$ iterations, will approximate the perimeter of the circle, measuring $2 \pi$. What this means is that $\lim_{n \to \infty} P_n = 2 \pi $.
Also note that every time we create a new polygon the number of sides doubles. Thus, after n steps we have a $2^n$ sided polygon and $P_n$ is:
$ P_n = 2^n \times s_n $
Manipulating the two expressions above we get:
$$ s_{n+1} = \sqrt{ 2 - \sqrt{4 - s_n^2} } $$Since we started with a square we have: $s_2 = \sqrt 2$. We can also consider $s_1 = 2$ representing a diameter line.
So, with this we have all equations needed to iteratively approximate $\pi$. Let's start by coding a function that gives us $s_{n+1}$:
In [1]:
from math import sqrt, pi
def side_next(side):
return sqrt(2. - sqrt(4. - side**2.0))
Our function aprox_pi() will compute the aproximation of $\pi$ for a $2^n$ polygon:
In [2]:
def aprox_pi(n = 10):
s = 2.0
for i in range(1, n):
s = side_next(s)
return 2.0**(n-1) * s
And the result is:
In [3]:
a_pi = aprox_pi()
a_pi
Out[3]:
So, in 10 iterations we got a very good approximation of $\pi$ using a 1024-sided polygon. (Note that since $P_\infty \rightarrow 2 \pi$ we need to divide the final result by two. aprox_pi() is automatically doing so.)
The error of the result is:
In [4]:
abs(pi - a_pi)
Out[4]:
That's Interesting. Let's see how good are the approximations generated by aprox_pi().
In [5]:
print("%10s \t %10s \t %20s \t %10s" % ("i", "Sides", "Pi", "Error"))
print("===================================================================")
for i in range(1, 31):
sides = 2.**i
a_pi = aprox_pi(i)
err = abs(pi - a_pi)
print("%10d \t %10d \t %20.10f \t %10.2e" % (i, sides, a_pi, err))
So, what's going on? We should expect better approximations as the number of sides increases, right? But, the best result we get is with 14 iterations and a polygon of 16384 sides. After that the approximations of $\pi$ get worse.
The problem is that our algorithm is not very good in terms of producing the end result. If you look at the expression $P_n = 2^n \times s_n$ what we are doing is multiplying a very large number (the number of sides) by a very small number (the length of a side). After a certain point, because we are using floating point, this is a recipe for disaster. In particular, for a 16384-sided polygon, the length of a side is approximatly:
In [6]:
2*pi/16384
Out[6]:
On a final note. Archimedes didn't have access to computers nor sofisticated ways of calculating square roots or sofisticated algebra. While in this notebook I've used his ideas of inscribing a polygon inside of a circle, his method was a little different in terms of calculation. He started with an hexagon (not a square) and ended-up with a 96-sided polygon. At the same time he ended up calculating his approximation of $\pi$ using fractions (which I still think is impressive):
"The ratio of the circumference of any circle to its diameter is greater than $3\tfrac{10}{71}$ but less than $3\tfrac{1}{7}$"
He did this by using rational approximations of certain square roots. Actually, no one knows where he got his roots! Check this page for a good description of this process of finding out $\pi$.
In [ ]: